An exploration of effective and useful visualisation of optical satellite images¶
Make your satellite image pop!¶
In [1]:
import re
import numpy as np
import matplotlib.pyplot as plt
import boto3
import rasterio
from rasterio.plot import show, show_hist
from skimage.exposure import equalize_adapthist
In [3]:
session = boto3.Session(profile_name='default')
In [4]:
scenes = [
'LC08_L1TP_205021_20180625_20180704_01_T1',
'LC08_L1TP_205021_20180727_20180731_01_T1'
]
In [5]:
def get_aws_band(scene, band, profile=False):
path, row = re.match(r'LC08_L1TP_(\d{3})(\d{3})', scene).groups()
prefix = f'c1/L8/{path}/{row}/{scene}'
with rasterio.Env(session=session):
with rasterio.open(
f's3://landsat-pds/{prefix}/{scene}_B{band}.TIF') as src:
if profile:
return src.profile
return src.read(out_shape=(src.height//10, src.width//10))
In [6]:
stacked_scenes = []
for scene in scenes:
red = get_aws_band(scene, 4)[0]
green = get_aws_band(scene, 3)[0]
blue = get_aws_band(scene, 2)[0]
stacked_scenes.append(np.stack([red, green, blue]))
Let's plot our downloaded image¶
We will use rasterio's plotting function (which wraps matplotlib)
In [11]:
fig, axes = plt.subplots(1,2, figsize=(30,15))
show(stacked_scenes[0], title="Scene 1", ax=axes[0])
show(stacked_scenes[1], title="Scene 2", ax=axes[1])
Out[11]:
In [12]:
fig, axes = plt.subplots(1,2, figsize=(30,10))
axes[0].set_xlim(0,65535)
axes[1].set_xlim(0,65535)
axes[0].set_ylim(0,140000)
axes[1].set_ylim(0,140000)
show_hist(
stacked_scenes[0], bins=50, lw=0.0, stacked=False, alpha=0.3,
histtype='stepfilled', title="Histogram", ax=axes[0])
show_hist(
stacked_scenes[1], bins=50, lw=0.0, stacked=False, alpha=0.3,
histtype='stepfilled', title="Histogram", ax=axes[1])
In [13]:
!gdalinfo data/LC08_L1TP_205021_20180625_20180704_01_T1/full_stack.tif
In [14]:
!gdalinfo data/LC08_L1TP_205021_20180625_20180704_01_T1/full_stack.tif | grep "Band"
Ok, let's just make it 8bit, that'll work, right...¶
Most plotting tools work with 8 bit images (we have 16 bit):
- Red; 256 tonal levels
- Green; 256 tonal levels
- Blue; 256 tonal levels
In [15]:
fig, axes = plt.subplots(1,2, figsize=(30,15))
show(stacked_scenes[0].astype(np.uint8), title="Scene 1", ax=axes[0])
show(stacked_scenes[1].astype(np.uint8), title="Scene 2", ax=axes[1])
Out[15]:
In [16]:
fig, axes = plt.subplots(1,2, figsize=(30,10))
axes[0].set_xlim(0,256)
axes[1].set_xlim(0,256)
axes[0].set_ylim(0,50000)
axes[1].set_ylim(0,50000)
show_hist(
stacked_scenes[0].astype(np.uint8), bins=50, lw=0.0, stacked=False,
alpha=0.3, histtype='stepfilled', title="Histogram", ax=axes[0])
show_hist(
stacked_scenes[1].astype(np.uint8), bins=50, lw=0.0, stacked=False,
alpha=0.3, histtype='stepfilled', title="Histogram", ax=axes[1])
That's pretty psychedelic, but not helpful¶
We have to be a bit smarter than that, we need to rescale!
In [17]:
def linear_rescale(image, in_range=[0, 65535], out_range=[1, 255]):
"""
Linear rescaling
"""
image = np.clip(image, in_range[0], in_range[1]) - in_range[0]
image = image / float(in_range[1] - in_range[0])
return (
image * (out_range[1] - out_range[0]) + out_range[0]
).astype(np.uint8)
In [18]:
linear_rescaled_images = [linear_rescale(stack) for stack in stacked_scenes]
In [19]:
fig, axes = plt.subplots(1,2, figsize=(30,15))
show(linear_rescaled_images[0], title="Scene 1", ax=axes[0])
show(linear_rescaled_images[1], title="Scene 2", ax=axes[1])
Out[19]:
In [20]:
fig, axes = plt.subplots(1,2, figsize=(30,10))
axes[0].set_xlim(0,256)
axes[1].set_xlim(0,256)
axes[0].set_ylim(0,120000)
axes[1].set_ylim(0,120000)
show_hist(
linear_rescaled_images[0], bins=50, lw=0.0, stacked=False, alpha=0.3,
histtype='stepfilled', title="Histogram", ax=axes[0])
show_hist(
linear_rescaled_images[1], bins=50, lw=0.0, stacked=False, alpha=0.3,
histtype='stepfilled', title="Histogram", ax=axes[1])
OK, Landsat 8 doesn't actually make use of the full 16 bits, so we will rescale within the range we have reasonable values:
In [21]:
linear_rescaled_images = [
linear_rescale(stack, in_range=[0, 14000]) for stack in stacked_scenes
]
In [22]:
fig, axes = plt.subplots(1,2, figsize=(30,15))
show(linear_rescaled_images[0], title="Scene 1", ax=axes[0])
show(linear_rescaled_images[1], title="Scene 2", ax=axes[1])
Out[22]:
In [23]:
fig, axes = plt.subplots(1,2, figsize=(30,10))
axes[0].set_xlim(0,256)
axes[1].set_xlim(0,256)
axes[0].set_ylim(0,120000)
axes[1].set_ylim(0,120000)
show_hist(
linear_rescaled_images[0], bins=50, lw=0.0, stacked=False, alpha=0.3,
histtype='stepfilled', title="Histogram", ax=axes[0])
show_hist(
linear_rescaled_images[1], bins=50, lw=0.0, stacked=False, alpha=0.3,
histtype='stepfilled', title="Histogram", ax=axes[1])
That's not very pretty¶
We can make better use of the colour range available, so lets stretch it.
In [24]:
def percentile_stretch(image, lower_perc=2, upper_perc=98):
mask = np.all(image != 0, axis=0).astype(np.uint8) * 255
image_out = np.zeros_like(image).astype(np.uint8)
for band in range(image.shape[0]):
input_range = np.percentile(
image[band][mask > 0],
(lower_perc, upper_perc)
).tolist()
image_out[band] = np.where(
mask,
linear_rescale(image[band],
in_range=input_range,
out_range=[0, 255]),
0
)
return image_out
In [25]:
percentile_stretched_images = [
percentile_stretch(stack) for stack in stacked_scenes
]
In [26]:
fig, axes = plt.subplots(1, 2, figsize=(30,15))
show(percentile_stretched_images[0], title="Scene 1", ax=axes[0])
show(percentile_stretched_images[1], title="Scene 2", ax=axes[1])
Out[26]:
In [27]:
fig, axes = plt.subplots(1,2, figsize=(30,10))
axes[0].set_xlim(0,256)
axes[1].set_xlim(0,256)
axes[0].set_ylim(0,80000)
axes[1].set_ylim(0,80000)
show_hist(
percentile_stretched_images[0], bins=50, lw=0.0, stacked=False, alpha=0.3,
histtype='stepfilled', title="Histogram", ax=axes[0])
show_hist(
percentile_stretched_images[1], bins=50, lw=0.0, stacked=False, alpha=0.3,
histtype='stepfilled', title="Histogram", ax=axes[1])
Give me more detail¶
That worked nicely with our clearish scene, but in the cloudy one, well the clouds look really interesting.
Contrast is lacking elsewhere...
Introducing Contrast Limited Adaptive histogram Equalization
In [28]:
def clahe_equalize(image, kernel_size=2000, clip_limit=0.03):
image_out = np.zeros_like(image).astype(np.float64)
for band in range(image.shape[0]):
image_out[band] = equalize_adapthist(
image[band], kernel_size=kernel_size, clip_limit=clip_limit
)
return image_out
In [29]:
raster_equalized_images = [
clahe_equalize(stack, clip_limit=0.05) for stack in stacked_scenes
]
In [30]:
fig, axes = plt.subplots(1,2, figsize=(30,15))
show(
(raster_equalized_images[0]* 255).astype(np.uint8),
title="Scene 1", ax=axes[0]
)
show(
(raster_equalized_images[1]* 255).astype(np.uint8),
title="Scene 2", ax=axes[1]
)
Out[30]:
In [31]:
fig, axes = plt.subplots(1,2, figsize=(30,10))
axes[0].set_xlim(0,256)
axes[1].set_xlim(0,256)
axes[0].set_ylim(0,80000)
axes[1].set_ylim(0,80000)
show_hist(
(raster_equalized_images[0]* 255).astype(np.uint8),
bins=50, lw=0.0, stacked=False, alpha=0.3, histtype='stepfilled',
title="Histogram", ax=axes[0]
)
show_hist(
(raster_equalized_images[1]* 255).astype(np.uint8),
bins=50, lw=0.0, stacked=False, alpha=0.3, histtype='stepfilled',
title="Histogram", ax=axes[1]
)
We could also remove the atmosphere...kind of...¶
The scenes still look very different, so not great for applying to all scenes, let's think about this in a physical sense.
We could try and remove the atmosphere, and ideally with an image based method (simpler, quicker, and no need for additional external data)
In [32]:
def find_dark_object_value(arr):
"""Find the value of the dark object
ie the first dark value that is not an outlier
"""
preval = None
step = arr.max() / 255.0
for val in np.unique(arr)[:100]:
if val == 0:
continue
if preval is not None and (val - preval) < step:
break
else:
preval = val
return preval
def dos(image):
image_out = np.zeros_like(image).astype(np.uint8)
# get the dark object for each band
for i, band in enumerate(image):
dark = find_dark_object_value(band)
image_out[i] = band - dark.astype(np.float)
image_out[image_out < 0] = 0
return image_out.astype(np.uint8)
In [33]:
prepped_scenes = [
linear_rescale(stack, [0, 14000]) for stack in stacked_scenes
]
dos_scenes = [dos(prepped_scene) for prepped_scene in prepped_scenes]
In [34]:
fig, axes = plt.subplots(1,2, figsize=(30,15))
show(linear_rescaled_images[0], title="Scene 1", ax=axes[0])
show(dos_scenes[0], title="Scene 1 - DOS", ax=axes[1])
Out[34]:
In [35]:
fig, axes = plt.subplots(1,2, figsize=(30,15))
show(linear_rescaled_images[1], title="Scene 2", ax=axes[0])
show(dos_scenes[1], title="Scene 2 - DOS", ax=axes[1])
Out[35]:
In [36]:
fig, axes = plt.subplots(1,2, figsize=(30,15))
show(dos_scenes[0], title="Scene 1 - DOS", ax=axes[0])
show(dos_scenes[1], title="Scene 2 - DOS", ax=axes[1])
Out[36]:
In [37]:
fig, axes = plt.subplots(1,2, figsize=(30,10))
axes[0].set_xlim(0,256)
axes[1].set_xlim(0,256)
axes[0].set_ylim(0,120000)
axes[1].set_ylim(0,120000)
show_hist(
dos_scenes[0], bins=50, lw=0.0, stacked=False, alpha=0.3,
histtype='stepfilled', title="Histogram", ax=axes[0])
show_hist(
dos_scenes[1], bins=50, lw=0.0, stacked=False, alpha=0.3,
histtype='stepfilled', title="Histogram", ax=axes[1])
Export for display¶
In [38]:
def write_rendered_rgb(rgb, profile, outname_prefix):
"""
Write a rendered RGB to GeoTIFF.
"""
profile.update(dtype=rasterio.uint8, count=3, photometric='RGB')
outname = f'data/{outname_prefix}_rendered.tif'
with rasterio.open(outname, 'w', **profile) as dst:
for k, arr in enumerate(rgb):
dst.write_band(k+1, arr)
In [39]:
profile = get_aws_band(scenes[1], 4, True)
equalized_image_8bit = (raster_equalized_images[1] * 255).astype(np.uint8)
write_rendered_rgb(equalized_image_8bit, profile, "scene2")
In [41]:
!gdalinfo data/scene2_rendered.tif | grep Band